Variation in vocal behavior under varying stress levels
Budgie INBRE stress and contact calls
Marcelo Araya-Salas, PhD & Timothy Wright, PhD
28-02-2022
Load packages
## add 'developer' to packages to be installed from github
x <- c("remotes", "lubridate", "readxl", "pbapply", "viridis", "ggplot2", "kableExtra", "knitr", "formatR", "randomForest", "Rtsne", "MASS", "adehabitatHR", "sp", "GGally", "spatstat", "raster", "vegan", "brms", "lme4", "dplyr", "purrr", "forcats", "tidyr", "modelr", "tidybayes", "cowplot", "ggrepel", "posterior", "ggridges", "maRce10/PhenotypeSpace")
aa <- lapply(x, function(y) {
# get pakage name
pkg <- strsplit(y, "/")[[1]]
pkg <- pkg[length(pkg)]
# check if installed, if not then install
if (!pkg %in% installed.packages()[,"Package"]) {
if (grepl("/", y)) remotes::install_github(y, force = TRUE) else
install.packages(y)
}
# load package
a <- try(require(pkg, character.only = T), silent = T)
if (!a) remove.packages(pkg)
})Functions and global parameters
opts_knit$set(root.dir = "..")
# set evaluation false
opts_chunk$set(fig.width = 10, fig.height = 6, warning = FALSE, message = FALSE, tidy = TRUE)
read_excel_df <- function(...) data.frame(read_excel(...))
trait_space_plot <- function(X, dimensions, group, indices, basecex = 1, title = "", colors = c("#3E4A89FF", "#35B779FF"), point.alpha = 0.7, background.indices = NULL, pch = 1, labels = c("sub-space", "total space"), legend = TRUE) {
total_coors <- as.ppp(as.matrix(X[, dimensions]), c(range(X[, dimensions[1]]), range(X[, dimensions[2]])))
total_space <- raster(density.ppp(total_coors))
xlim <- range(X[, dimensions[1]])
ylim <- range(X[, dimensions[2]])
# ylim <- if (is.null(background.indices))
# range(X[, dimensions[2]]) else
# c(min(X[, dimensions[2]]), c(max(X[, dimensions[2]]) + (max(X[, dimensions[2]]) - min(X[, dimensions[2]])) * 0.11))
plot(x = X[, dimensions[1]], y = X[, dimensions[2]], col = "white", cex = basecex, xlab = dimensions[1], ylab= dimensions[2], xlim = xlim, ylim = ylim, cex.lab = basecex, xaxs="i", yaxs="i")
# add background group
if (!is.null(background.indices)){
# get points
Y_bg <- X[background.indices, ]
if (nrow(Y_bg) > 1 & point.alpha > 0)
points(Y_bg[, dimensions[1]], Y_bg[, dimensions[2]], col = adjustcolor(colors[2], alpha.f = point.alpha), pch = pch)
if (nrow(Y_bg) > 4){
coors <- as.ppp(as.matrix(Y_bg[, dimensions]), c(range(Y_bg[, dimensions[1]]), range(Y_bg[, dimensions[2]])))
raster_dens <- raster(density.ppp(coors))
palpre <- colorRampPalette(c("white", colors[2]))
colsn <- 10
palpre2 <- sapply(1:colsn, function(x) adjustcolor(col = palpre(colsn)[x], alpha = seq(0.1, 0.9, length.out = 20)[x]))
image(raster_dens, add = TRUE, col = palpre2)
}
}
Y <- X[indices, ]
if (nrow(Y) > 1 & point.alpha > 0)
points(Y[, dimensions[1]], Y[, dimensions[2]], col = adjustcolor(colors[1], alpha.f = point.alpha), pch = pch)
if (nrow(Y) > 4){
coors <- as.ppp(as.matrix(Y[, dimensions]), c(range(Y[, dimensions[1]]), range(Y[, dimensions[2]])))
raster_dens <- raster(density.ppp(coors))
palpre <- colorRampPalette(c("white", colors[1]))
colsn <- 10
palpre2 <- sapply(1: colsn, function(x) adjustcolor(col = palpre(colsn)[x], alpha = seq(0.1, 0.9, length.out = 20)[x]))
image(raster_dens, add = TRUE, col = palpre2)
}
usr <- par()$usr
# legend
if (!is.null(background.indices) & legend)
# legend(y = rep(usr[4] + ((usr[4] - usr[3]) * 0.05), 2), x = c(usr[1] + (usr[2] - usr[1]) * 0.25, usr[1] + (usr[2] - usr[1]) * 0.75), col = colors, legend = labels, pch = c(pch, pch), cex = basecex * 1.5, bty = "n", horiz = TRUE, xjust = 0)
legend("topright", col = colors, legend = labels, pch = c(pch, pch), cex = basecex * 1.5, xjust = 0, bg = adjustcolor("white", alpha.f = 0.5))
title(title, cex.main = basecex * 1.2)
}
# for reading months in english format
sl <- Sys.setlocale(locale = "en_US.UTF-8")
cols <- viridis(10, alpha = 0.7)
summary_brm_model <- function(x){
summ <- summary(x)$fixed
fit <- x$fit
betas <- grep("^b_", names(fit@sim$samples[[1]]), value = TRUE)
hdis <- t(sapply(betas, function(y) hdi(fit@sim$samples[[1]][[y]]))
)
summ_table <- data.frame(summ, hdis, iterations = attr(fit, "stan_args")[[1]]$iter, chains = length(attr(fit, "stan_args")))
summ_table <- summ_table[rownames(summ_table) != "Intercept", c("Estimate", "Rhat", "Bulk_ESS", "l.95..CI", "u.95..CI", "iterations", "chains")]
summ_table <- as.data.frame(summ_table)
summ_table$Rhat <- round(summ_table$Rhat, digits = 3)
summ_table$CI_low <- round(unlist(summ_table$l.95..CI), digits = 3)
summ_table$CI_high <- round(unlist(summ_table$u.95..CI), digits = 3)
summ_table$l.95..CI <- summ_table$u.95..CI <- NULL
out <- lapply(betas, function(y) data.frame(variable = y, value = sort(fit@sim$samples[[1]][[y]], decreasing = FALSE)))
posteriors <- do.call(rbind, out)
posteriors <- posteriors[posteriors$variable != "b_Intercept", ]
gg <- ggplot(data = posteriors, aes(y = variable, x = value, fill = stat(quantile))) +
geom_density_ridges_gradient(quantile_lines = TRUE, quantile_fun = HDInterval::hdi, vline_linetype = 2) +
theme_classic() +
xlim(c(min(summ_table[ , "CI_low"]), max(summ_table[ , "CI_high"]))) +
geom_vline(xintercept = 0, col = "red") +
scale_fill_manual(values = c("transparent", "lightblue", "transparent"), guide = "none") + ggtitle(x$formula)
summ_table$Rhat <- ifelse(summ_table$Rhat > 1.05, cell_spec(summ_table$Rhat, "html", color ="white", background = "red", bold = TRUE, font_size = 12), cell_spec(summ_table$Rhat, "html"))
signif <- summ_table[,"CI_low"] * summ_table[,"CI_high"] > 0
df1 <- kbl(summ_table, row.names = TRUE, escape = FALSE, format = "html", digits = 3)
df1 <- row_spec(kable_input = df1,row = which(summ_table$CI_low * summ_table$CI_high > 0), background = adjustcolor(cols[9], alpha.f = 0.3))
df1 <- kable_styling(df1, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 12)
print(df1)
print(gg)
}# read ext sel tab calls
sels <- read.csv("./data/processed/tailored_budgie_calls_sel_tab.csv")
# keep only spectrographic parameters
sels <- sels[, c("sound.files", "selec", "duration", "meanfreq", "sd", "freq.median",
"freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom",
"mindom", "maxdom", "dfrange", "modindx", "meanpeakf")]
sels$ID <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][1])
sels$month <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][2])
sels$day <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][3])
sels$year <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][4])
sels$date <- paste(sels$day, substr(sels$month, 0, 3), sels$year, sep = "-")
sels$date <- as.Date(sels$date, format = "%d-%b-%Y")Counts per individual
df <- as.data.frame(table(sels$ID))
names(df) <- c("ID", "Sample_size")
df <- df[order(df$Sample_size, decreasing = FALSE), ]
kb <- kable(df, row.names = FALSE)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)| ID | Sample_size |
|---|---|
| 132YGMM | 6 |
| 125YGMM | 12 |
| 160YGHM | 16 |
| 310YGHM | 21 |
| 393YGHM | 25 |
| 303YGHM | 37 |
| 398WBHM | 41 |
| 365YGHM | 44 |
| 399YGLM | 46 |
| 300YGHM | 47 |
| 270YGHM | 64 |
| 407YGHM | 97 |
| 386WBHM | 100 |
| 377WWLM | 108 |
| 367WWNM | 119 |
| 371YYLM | 149 |
| 397YGHM | 175 |
| 378YYLM | 195 |
| 404WBHM | 196 |
| 258YGHM | 223 |
| 389WWLM | 230 |
| 262YGHM | 306 |
| 400YGHM | 360 |
| 388YGLM | 377 |
| 327YYLM | 404 |
| 396YBHM | 455 |
| 403WBHM | 566 |
| 356WBHM | 761 |
| 361YGHM | 767 |
| 402YGHM | 849 |
| 395WBHM | 854 |
| 360YGHM | 900 |
| 390WBHM | 935 |
| 405YBHM | 975 |
| 385YBHM | 1256 |
| 394YBHM | 1693 |
Add metadata
metadat <- read_excel_df("./data/raw/INBREStress_MasterDataSheet_14Nov19.xlsx")
# head(metadat)
sels$ID[sels$ID == "125YGMM"] <- "125YGHM"
sels$ID[sels$ID == "394YBHM"] <- "394WBHM"
# setdiff(sels$ID, metadat$Bird.ID) setdiff(metadat$Bird.ID, sels$ID)
sels$treatment <- sapply(1:nrow(sels), function(x) {
metadat$Treatment[metadat$Bird.ID == sels$ID[x]][1]
})
sels$treatment.room <- sapply(1:nrow(sels), function(x) {
metadat$Treatment.Room[metadat$Bird.ID == sels$ID[x]][1]
})
sels$round <- sapply(1:nrow(sels), function(x) {
metadat$Round[metadat$Bird.ID == sels$ID[x]][1]
})
sels$source.room <- sapply(1:nrow(sels), function(x) {
metadat$Source.Room[metadat$Bird.ID == sels$ID[x]][1]
})
sels$record.group <- sapply(1:nrow(sels), function(x) {
metadat$Record.Group[metadat$Bird.ID == sels$ID[x]][1]
})
# add week
out <- lapply(unique(sels$round), function(x) {
Y <- sels[sels$round == x, ]
min_date <- min(Y$date)
week_limits <- min_date + seq(0, 100, by = 7)
Y$week <- NA
for (i in 2:length(week_limits)) Y$week[Y$date >= week_limits[i - 1] & Y$date <
week_limits[i]] <- i - 1
return(Y)
})
sels <- do.call(rbind, out)
sels$cort.baseline <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$cort.stress <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$stress.response <- sels$cort.stress - sels$cort.baseline
sels$weight <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$breath.count <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
# remove week 5
sels <- sels[sels$week != 5, ]Exploratory graph
Physiological parameters
breath.count <- stack(metadat[, c("D3.Breath.Count", "D7.Breath.Count", "D14.Breath.Count",
"D21.Breath.Count", "D28.Breath.Count")])
weight <- stack(metadat[, c("D3.Bird.Weight..g.", "D7.Bird.Weight..g.", "D14.Bird.Weight..g.",
"D21.Bird.Weight..g.", "D28.Bird.Weight..g.")])
cort.stress <- stack(metadat[, c("D3.CORT.Stress", "D7.CORT.Stress", "D14.CORT.Stress",
"D21.CORT.Stress", "D28.CORT.Stress")])
cort.baseline <- stack(metadat[, c("D3.CORT.Baseline", "D7.CORT.Baseline", "D14.CORT.Baseline",
"D21.CORT.Baseline", "D28.CORT.Baseline")])
stress <- data.frame(metadat[, c("Bird.ID", "Treatment", "Round", "Treatment.Room")],
week = breath.count$ind, breath.count = breath.count$values, weight = weight$values,
cort.stress = cort.stress$values, cort.baseline = cort.baseline$values, stress.response = cort.stress$values -
cort.baseline$values)
# head(stress)
stress$week <- factor(sapply(strsplit(as.character(stress$week), "\\."), "[[", 1),
levels = c("D3", "D7", "D14", "D21", "D28"))
stress$day <- as.numeric(gsub("D", "", as.character(stress$week)))
stress$round <- paste("Round", stress$Round)
stress$treatment <- factor(stress$Treatment, levels = c("Control", "Medium Stress",
"High Stress"))
# remove 5th week
stress <- stress[stress$week != "D28", ]
stress_l <- lapply(stress$Bird.ID, function(x) {
X <- stress[stress$Bird.ID == x, ]
X$breath.count <- X$breath.count - X$breath.count[X$week == "D3"]
X$weight <- X$weight - X$weight[X$week == "D3"]
X$cort.stress <- X$cort.stress - X$cort.stress[X$week == "D3"]
X$cort.baseline <- X$cort.baseline - X$cort.baseline[X$week == "D3"]
X$stress.response <- X$stress.response - X$stress.response[X$week == "D3"]
return(X)
})
stress <- do.call(rbind, stress_l)
ggplot(data = stress, aes(x = day, y = weight, group = Bird.ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Weight (g)", x = "Day sample was taken") +
theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))
ggplot(data = stress, aes(x = day, y = breath.count, group = Bird.ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Breath count", x = "Day sample was taken") +
theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))
ggplot(data = stress, aes(x = day, y = cort.baseline, group = Bird.ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Baseline CORT", x = "Day sample was taken") +
theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))
ggplot(data = stress, aes(x = day, y = cort.stress, group = Bird.ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Stress CORT", x = "Day sample was taken") +
theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))
ggplot(data = stress, aes(x = day, y = stress.response, group = Bird.ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Stress response", x = "Day sample was taken") +
theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))
# stress_l <- lapply(unique(stress$Bird.ID), function(x){ X <-
# stress[stress$Bird.ID == x, ] X$cort.baseline.change <- X$cort.baseline - X$
# })
standard_error <- function(x) sd(x)/sqrt(length(x))
agg_stress <- aggregate(cort.baseline ~ week + treatment, stress, mean)
agg_stress$se <- aggregate(cort.baseline ~ week + treatment, stress, standard_error)[,
3]
agg_stress$Week <- 1:4
ggplot(data = agg_stress, aes(x = Week, y = cort.baseline, fill = treatment)) + geom_bar(stat = "identity") +
geom_errorbar(aes(ymin = cort.baseline - se, ymax = cort.baseline + se), width = 0.1) +
# geom_line(size = 1.2) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") +
labs(y = "Baseline CORT", x = "Week") + theme_classic(base_size = 30) + theme(legend.position = "none")
FOXP2
foxp2 <- read.csv("./data/raw/stress_IHC_counts.csv")
foxp2$Treatment[grep("HIGH", foxp2$Treatment)] <- "High stress"
foxp2$Treatment[grep("Control", foxp2$Treatment, ignore.case = T)] <- "Control"
foxp2$Treatment[grep("MED", foxp2$Treatment)] <- "Medium stress"
foxp2$Treatment <- factor(foxp2$Treatment, levels = c("Control", "Medium stress",
"High stress"))
ggplot(data = foxp2, aes(x = Treatment, y = FoxP2.Counts.MMSt.Striatum, fill = Treatment)) +
geom_boxplot() + geom_point() + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
labs(y = "% FoxP2 expression", x = "Treatment") + theme_classic(base_size = 30) +
# scale_x_discrete( labels = c(1:4)) +
theme(legend.position = "none")Behavioral parameters
Call counts per treatment
call_count <- aggregate(selec ~ week + treatment + round, data = sels, length)
call_count$Week <- factor(call_count$week, levels = (unique(call_count$week))[5:1])
call_count$round <- paste("Round", call_count$round)
call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
call_count$stress.response <- aggregate(selec ~ week + treatment + round, data = sels,
mean)[, 4]
ggplot(data = call_count, aes(x = treatment, y = selec, fill = Week)) + geom_bar(stat = "identity") +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + coord_flip() + facet_wrap(~round,
ncol = 2) + labs(y = "Number of calls", x = "Treatment") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.15))
call_count2 <- aggregate(selec ~ week + treatment + round + ID, data = sels, length)
agg_count <- aggregate(selec ~ week + treatment, call_count2, mean)
agg_count$se <- aggregate(selec ~ week + treatment, call_count2, standard_error)[,
3]
agg_count$treatment <- factor(agg_count$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
ggplot(data = agg_count, aes(x = week, y = selec, fill = treatment)) + geom_bar(stat = "identity") +
geom_errorbar(aes(ymin = selec - se, ymax = selec + se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Number of calls",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = "none")
Call counts per bird and treatment
# Call counts per treatment by individual (in color)
call_count <- aggregate(selec ~ week + treatment + round + ID, data = sels, length)
call_count$Week <- factor(call_count$week, levels = (1:5))
call_count$round <- paste("Round", call_count$round)
call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
ggplot(data = call_count, aes(x = week, y = selec, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Call count", x = "Week") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.25))
Acoustic parameters
call_week_params <- aggregate(cbind(duration, meanfreq, sd, freq.median, freq.IQR,
time.IQR, skew, kurt, sp.ent, time.ent, entropy, meandom, mindom, maxdom, dfrange,
modindx, meanpeakf) ~ week + treatment + round, data = sels, mean)
call_week_params_sd <- aggregate(cbind(duration, meanfreq, sd, freq.median, freq.IQR,
time.IQR, skew, kurt, sp.ent, time.ent, entropy, meandom, mindom, maxdom, dfrange,
modindx, meanpeakf) ~ week + treatment + round, data = sels, sd)
names(call_week_params_sd) <- names(call_week_params) <- c("week", "treatment", "round",
"Duration", "Mean_frequency_(kHz)", "SD_frequency_(kHz)", "Median_frequency_(kHz)",
"Interquantile_frequency_range_(kHz)", "Interquantile_time_range_(s)", "Skewness",
"Kurtosis", "Spectral_entropy", "Time_entropy", "Overall entropy", "Mean_dominant_frequency_(kHz)",
"Minimum_dominant_frequency_(kHz)", "Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)",
"Modulation_index", "Mean_peak_frequency_(kHz)")
call_week_params_sd <- call_week_params_sd[, c("Duration", "Mean_frequency_(kHz)",
"SD_frequency_(kHz)", "Median_frequency_(kHz)", "Interquantile_frequency_range_(kHz)",
"Interquantile_time_range_(s)", "Skewness", "Kurtosis", "Spectral_entropy", "Time_entropy",
"Overall entropy", "Mean_dominant_frequency_(kHz)", "Minimum_dominant_frequency_(kHz)",
"Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)", "Modulation_index",
"Mean_peak_frequency_(kHz)")]
names(call_week_params_sd) <- paste(names(call_week_params_sd), "sd", sep = ".")
call_week_params <- cbind(call_week_params, call_week_params_sd)
call_week_params$round <- paste("Round", call_week_params$round)
call_week_params$treatment <- factor(call_week_params$treatment, levels = c("Control",
"Medium Stress", "High Stress"))
pd <- position_dodge(0.3)
ggs <- lapply(c("Duration", "Mean_frequency_(kHz)", "SD_frequency_(kHz)", "Median_frequency_(kHz)",
"Interquantile_frequency_range_(kHz)", "Interquantile_time_range_(s)", "Skewness",
"Kurtosis", "Spectral_entropy", "Time_entropy", "Overall entropy", "Mean_dominant_frequency_(kHz)",
"Minimum_dominant_frequency_(kHz)", "Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)",
"Modulation_index", "Mean_peak_frequency_(kHz)"), function(i) {
call_week_params$var <- call_week_params[, i]
call_week_params$var.min <- call_week_params$var - call_week_params[, paste(i,
"sd", sep = ".")]
call_week_params$var.max <- call_week_params$var + call_week_params[, paste(i,
"sd", sep = ".")]
ggplot(call_week_params, aes(x = week, y = var, col = treatment)) + geom_point(size = 4,
position = pd) + geom_errorbar(aes(ymin = var.min, ymax = var.max), width = 0.3,
position = pd, size = 1.7) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = i, x = "Week") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.25))
})
ggs## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
##
## [[10]]
##
## [[11]]
##
## [[12]]
##
## [[13]]
##
## [[14]]
##
## [[15]]
##
## [[16]]
##
## [[17]]
Acoustic space projection
scale_acous_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median",
"freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom",
"mindom", "maxdom", "dfrange", "modindx", "meanpeakf")])
urfmod <- randomForest(x = scale_acous_param, importance = TRUE, proximity = TRUE,
ntree = 10000)
saveRDS(urfmod, "./data/processed/unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")
mds_urf_prox <- cmdscale(1 - urfmod$proximity, k = 2)
saveRDS(mds_urf_prox, "./data/processed/MDS_unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")Random forest
urfmod <- readRDS("./data/processed/unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")
mds_urf_prox <- readRDS("./data/processed/MDS_unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")
sels$MDS1 <- mds_urf_prox[, 1]
sels$MDS2 <- mds_urf_prox[, 2]
print(1 - sum(urfmod$votes[, 1] > urfmod$votes[, 2])/nrow(urfmod$votes))## [1] 0.005058366Proportion of real data classified as synthetic data by the unsupervised random forest: 0.005058
ggplot(sels, aes(x = MDS1, y = MDS2, col = as.factor(round))) + geom_point() + labs(color = "Round") +
scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 30) + theme(legend.position = c(0.62,
0.3)) + guides(colour = guide_legend(override.aes = list(size = 10)))
PCA
pca <- prcomp(x = sels[, c("duration", "meanfreq", "sd", "freq.median", "freq.IQR",
"time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom", "mindom",
"maxdom", "dfrange", "modindx", "meanpeakf")], scale. = TRUE)
Y <- as.data.frame(pca$x[, 1:2])
names(Y) <- c("PC1", "PC2")
sels <- data.frame(sels, Y)
ggplot(sels, aes(x = PC1, y = PC2, col = as.factor(round))) + geom_point() + labs(color = "Round") +
scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) + theme(legend.position = c(0.9,
0.8)) + guides(colour = guide_legend(override.aes = list(size = 10)))
t-SNE
scale_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median", "freq.IQR",
"time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom", "mindom",
"maxdom", "dfrange", "modindx", "meanpeakf")])
tsne <- Rtsne(scale_param, dims = 2, perplexity = 30, verbose = FALSE, max_iter = 5000)
saveRDS(tsne, "./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")tsne <- readRDS("./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")
Y <- as.data.frame(tsne$Y)
names(Y) <- c("TSNE1", "TSNE2")
sels <- data.frame(sels, Y)
ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(round))) + geom_point() +
labs(color = "Round") + scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) +
theme(legend.position = c(0.955, 0.25)) + guides(colour = guide_legend(override.aes = list(size = 10)))
sels$treatment <- factor(sels$treatment, levels = c("Control", "Medium Stress", "High Stress"))
ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(treatment))) + geom_point() +
labs(color = "Treatment") + scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) +
theme(legend.position = c(0.88, 0.2)) + guides(colour = guide_legend(override.aes = list(size = 10)))
Acoustic space by individual (each point) for the 3 dimension reduction methods
- Only individuals with at least 20 calls
- Using rarefaction with n = the minimum sample size across all individuals
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab[tab > 20]), ]
pca_areas <- rarefact_space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID",
parallel = 4, pb = TRUE, type = "density")
rf_mds_areas <- rarefact_space_size(X = Y, dimensions = c("MDS1", "MDS2"), group = "ID",
parallel = 4, pb = TRUE, type = "density")
tsne_areas <- rarefact_space_size(X = Y, dimensions = c("TSNE1", "TSNE2"), group = "ID",
parallel = 4, pb = TRUE, type = "density")
areas <- data.frame(pca_areas[, -ncol(pca_areas)], pca.area = pca_areas$mean.size,
rf.mds.area = rf_mds_areas$mean.size, tsne.area = tsne_areas$mean.size)
saveRDS(areas, "./data/processed/acoustic_space_area_by_individual_for_each_dimension_reduction_method.RDS")areas <- readRDS("./data/processed/acoustic_space_area_by_individual_for_each_dimension_reduction_method.RDS")
ggpairs(areas, columns = which(names(areas) %in% c("pca.area", "rf.mds.area", "tsne.area")),
columnLabels = c("PCA", "Rand.for. MDS", "t-SNE")) + theme_minimal(base_size = 30)
Individual acoustic space across weeks
- PCA and t-SNE used for dimensionality reduction
- Only weeks with at least 6 calls were included
PCA
sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID.week)
# sort(tab)
Y <- sels[sels$ID.week %in% names(tab)[tab >= 6], ]
# run with rarefaction
areas_by_week_raref <- rarefact_space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID.week",
parallel = 1, pb = TRUE, type = "density")
# run without rarefaction
areas_by_week <- space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID.week",
parallel = 1, pb = TRUE, type = "density", extent = 2)
areas_by_week$raref.area <- areas_by_week_raref$mean.size
# add metadata
areas_by_week$ID <- sapply(strsplit(areas_by_week$group, "-"), "[[", 1)
areas_by_week$week <- factor(sapply(strsplit(areas_by_week$group, "-"), "[[", 2),
levels = 1:5)
areas_by_week$treatment <- sapply(1:nrow(areas_by_week), function(x) metadat$Treatment[metadat$Bird.ID ==
areas_by_week$ID[x]][1])
areas_by_week$round <- sapply(1:nrow(areas_by_week), function(x) metadat$Round[metadat$Bird.ID ==
areas_by_week$ID[x]][1])
areas_by_week$round <- paste("Round", areas_by_week$round)
call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
areas_by_week$treatment <- factor(areas_by_week$treatment, levels = c("Control",
"Medium Stress", "High Stress"))
saveRDS(areas_by_week, "./data/processed/acoustic_space_area_by_individual_and_week.RDS")With rarefaction
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")
ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))
With rarefaction and compared to week 1
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")
areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {
X <- areas_by_week[areas_by_week$ID == x, ]
X$raref.area <- X$raref.area - X$raref.area[X$week == min(as.numeric(as.character(X$week)))]
return(X)
}))
ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))
Without rarefaction
ggplot(data = areas_by_week[areas_by_week$week != 5, ], aes(x = week, y = size, group = ID,
col = treatment)) + geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))
Without rarefaction and compared to week 1
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")
areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {
X <- areas_by_week[areas_by_week$ID == x, ]
X$size <- X$size - X$size[X$week == min(as.numeric(as.character(X$week)))]
return(X)
}))
ggplot(data = areas_by_week, aes(x = week, y = size, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))
Acoustic space area by number of calls per individual and week
ggplot(data = areas_by_week, aes(x = n, y = size, col = treatment)) + geom_point() +
# geom_smooth(size = 1.2, method = 'lm') +
scale_color_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~round, ncol = 2, scales = "free_x") +
labs(y = "Acoustic space area", x = "Number of calls (n)") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.25))
- Acoustic space area by number of calls per individual and week by treatment
- Must have at least 2 calls
ggplot(data = areas_by_week, aes(x = n, y = log(size + 10), col = treatment)) + geom_point() +
geom_smooth(size = 1.2, method = "lm") + scale_color_viridis_d(begin = 0.1, end = 0.9) +
# facet_wrap(~ round, ncol = 2, scales = 'free_x') +
labs(y = "log(acoustic space area)", x = "Number of calls (n)") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.7))
Save acoustic space plots
trait_space_plot(X = sels, indices = which(sels$ID == "395WBHM"), dimensions = c("TSNE1",
"TSNE2"), point.alpha = 0.2, pch = 1)
#
trait_space_plot(X = sels, indices = which(sels$ID == "395WBHM"), dimensions = c("TSNE1",
"TSNE2"), background.indices = which(sels$ID != "395WBHM" & sels$record.group ==
sels$record.group[sels$ID == "395WBHM"][1]), point.alpha = 0.2, pch = 1)# Y <- sels[sels$ID == '395WBHM',] Y <- sels[sels$ID %in% c('395WBHM',
# '394WBHM', '360YGHM', '385YBHM'), ]
# trait_overlap(Y, dimensions = c('TSNE1', 'TSNE2'), group = 'ID', type =
# 'density', overlap.type = 'assymetric')
# X = Y dimensions = c('TSNE1', 'TSNE2') level <- '395WBHM' basecex <- 1 group
# <- 'ID' indices <- which(X[, group] == level)
trait_space_plot(X = sels, indices = which(sels$ID == "395WBHM"), dimensions = c("TSNE1",
"TSNE2"), background.indices = which(sels$ID != "395WBHM" & sels$record.group ==
sels$record.group[sels$ID == "395WBHM"][1]), point.alpha = 0.2, pch = 1)
# with points
out <- pblapply(unique(sels$ID), function(x) {
tiff(filename = file.path("./images", paste(x, sels$treatment[sels$ID == x][1],
"v2.tiff", sep = "-")), res = 120, width = 800, height = 800)
par(mfrow = c(2, 2))
# W <- sels[sels$ID == x, ]
for (i in 1:4) try(trait_space_plot(X = sels, dimensions = c("TSNE1", "TSNE2"),
indices = which(sels$ID == x & sels$week == i), background.indices = which(sels$ID !=
x & sels$week == i & sels$record.group == sels$record.group[sels$ID ==
x][1]), title = paste(x, "week", i), basecex = 1, labels = c(x, "group"),
legend = FALSE), silent = TRUE)
dev.off()
})TSNE
sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID.week)
# sort(tab)
Y <- sels[sels$ID.week %in% names(tab)[tab >= 6], ]
# run with rarefaction
areas_by_week_raref <- rarefact_space_size(X = Y, dimensions = c("TSNE1", "TSNE2"),
group = "ID.week", parallel = 1, pb = TRUE, type = "density")
# run without rarefaction
areas_by_week <- space_size(X = Y, dimensions = c("TSNE1", "TSNE2"), group = "ID.week",
parallel = 4, pb = TRUE, type = "density")
areas_by_week$raref.area <- areas_by_week_raref$mean.size
# add metadata
areas_by_week$ID <- sapply(strsplit(areas_by_week$group, "-"), "[[", 1)
areas_by_week$week <- factor(sapply(strsplit(areas_by_week$group, "-"), "[[", 2),
levels = 1:5)
areas_by_week$treatment <- sapply(1:nrow(areas_by_week), function(x) metadat$Treatment[metadat$Bird.ID ==
areas_by_week$ID[x]][1])
areas_by_week$round <- sapply(1:nrow(areas_by_week), function(x) metadat$Round[metadat$Bird.ID ==
areas_by_week$ID[x]][1])
areas_by_week$round <- paste("Round", areas_by_week$round)
call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
areas_by_week$treatment <- factor(areas_by_week$treatment, levels = c("Control",
"Medium Stress", "High Stress"))
saveRDS(areas_by_week, "./data/processed/acoustic_space_area_by_individual_and_week_TSNE.RDS")With rarefaction
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week_TSNE.RDS")
ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))
With rarefaction and compared to week 1
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week_TSNE.RDS")
areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {
X <- areas_by_week[areas_by_week$ID == x, ]
X$raref.area <- X$raref.area - X$raref.area[X$week == min(as.numeric(as.character(X$week)))]
return(X)
}))
ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))
Without rarefaction
ggplot(data = areas_by_week[areas_by_week$week != 5, ], aes(x = week, y = size, group = ID,
col = treatment)) + geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))
Without rarefaction and compared to week 1
areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {
X <- areas_by_week[areas_by_week$ID == x, ]
X$size <- X$size - X$size[X$week == min(as.numeric(as.character(X$week)))]
return(X)
}))
ggplot(data = areas_by_week, aes(x = week, y = size, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))
Acoustic space area by number of calls per individual and week
ggplot(data = areas_by_week, aes(x = n, y = size, col = treatment)) + geom_point() +
# geom_smooth(size = 1.2, method = 'lm') +
scale_color_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~round, ncol = 2, scales = "free_x") +
labs(y = "Acoustic space area", x = "Number of calls (n)") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.25))
- Acoustic space area by number of calls per individual and week by treatment
- Must have at least 2 calls
ggplot(data = areas_by_week, aes(x = n, y = log(size + 10), col = treatment)) + geom_point() +
geom_smooth(size = 1.2, method = "lm") + scale_color_viridis_d(begin = 0.1, end = 0.9) +
# facet_wrap(~ round, ncol = 2, scales = 'free_x') +
labs(y = "log(acoustic space area)", x = "Number of calls (n)") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.7))
Acoustic space overlap
- 2 methods: pairwise distance between observations and kernel density overlap
Individual overlap to its own group
PCA
sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 2], ]
group_ovlp_l <- lapply(unique(Y$ID.week), function(x) {
# group data
X <- Y[Y$record.group == Y$record.group[Y$ID.week == x][1] & Y$week == Y$week[Y$ID.week ==
x][1], ]
# label all other individuals as group
X$ID2 <- ifelse(X$ID.week == x, "focal", "group")
# run with density
ovlp <- if (min(table(X$ID2)) > 4)
rarefact_space_similarity(X = X, dimensions = c("PC1", "PC2"), group = "ID2",
parallel = 1, pb = FALSE, outliers = 0.95, type = "mean.density.overlap") else matrix(rep(NA, 3), nrow = 1)
dsts <- if (min(table(X$ID2)) > 1)
rarefact_space_similarity(X = X, dimensions = c("PC1", "PC2"), group = "ID2",
parallel = 1, pb = FALSE, outliers = 0.95, type = "distance") else matrix(rep(NA, 3), nrow = 1)
out <- data.frame(ID = Y$ID[Y$ID.week == x][1], group = Y$record.group[Y$ID.week ==
x][1], week = Y$week[Y$ID.week == x][1], overlap.to.group = ovlp[1, 3], distance.to.group = dsts[1,
3], treatment = Y$treatment[Y$ID.week == x][1], round = Y$round[Y$ID.week ==
x][1])
return(out)
})
group_ovlp <- do.call(rbind, group_ovlp_l)
saveRDS(group_ovlp, "./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")
group_ovlp <- do.call(rbind, lapply(unique(group_ovlp$ID), function(x) {
X <- group_ovlp[group_ovlp$ID == x, ]
X$overlap.to.group <- X$overlap.to.group - X$overlap.to.group[X$week == min(as.numeric(as.character(X$week)))]
X$distance.to.group <- X$distance.to.group - X$distance.to.group[X$week == min(as.numeric(as.character(X$week)))]
return(X)
}))
group_ovlp$treatment <- factor(group_ovlp$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
ggplot(data = group_ovlp, aes(x = week, y = overlap.to.group, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to group acoustic space",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
0.25))
ggplot(data = group_ovlp, aes(x = week, y = overlap.to.group, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Individual distance to group acoustic space",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
0.25))
sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID.week)
Y <- sels[sels$ID.week %in% names(tab)[tab >= 6], ]
# run with rarefaction
overlap_by_week <- rarefact_space_similarity(X = Y, dimensions = c("PC1", "PC2"),
group = "ID.week", parallel = 1, pb = TRUE, outliers = 0.95, type = "mean.density.overlap")
saveRDS(overlap_by_week, "./data/processed/acoustic_space_overlap_by_individual_and_week.RDS")TSNE
sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 6], ]
group_ovlp_l <- lapply(unique(Y$ID.week), function(x) {
# group data
X <- Y[Y$record.group == Y$record.group[Y$ID.week == x][1] & Y$week == Y$week[Y$ID.week ==
x][1], ]
# label all other individuals as group
X$ID2 <- ifelse(X$ID.week == x, "focal", "group")
# run with density
ovlp <- if (min(table(X$ID2)) >= 6)
rarefact_space_similarity(X = X, dimensions = c("TSNE1", "TSNE2"), group = "ID2",
parallel = 1, pb = FALSE, type = "mean.density.overlap")[, 1:3] else matrix(rep(NA, 3), nrow = 1)
dsts <- if (min(table(X$ID2)) >= 2)
rarefact_space_similarity(X = X, dimensions = c("TSNE1", "TSNE2"), group = "ID2",
parallel = 1, pb = FALSE, type = "distance")[, 1:3] else matrix(rep(NA, 3), nrow = 1)
out <- data.frame(ID = X$ID[X$ID.week == x][1], group = X$record.group[X$ID.week ==
x][1], week = X$week[X$ID.week == x][1], overlap.to.group = ovlp[1, 3], distance.to.group = dsts[1,
3], treatment = X$treatment[X$ID.week == x][1], round = X$round[X$ID.week ==
x][1])
return(out)
})
group_ovlp <- do.call(rbind, group_ovlp_l)
saveRDS(group_ovlp, "./data/processed/acoustic_space_density_overlap_to_group_by_week_TSNE.RDS")group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week_TSNE.RDS")
# comparing to week 1 group_ovlp <- do.call(rbind,
# lapply(unique(group_ovlp$ID), function(x){ X <- group_ovlp[group_ovlp$ID ==
# x, ] X$overlap.to.group <- X$overlap.to.group / X$overlap.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] X$distance.to.group <-
# X$distance.to.group - X$distance.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] return(X) }))
group_ovlp$treatment <- factor(group_ovlp$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
ggplot(data = group_ovlp, aes(x = week, y = overlap.to.group, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to group acoustic space",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
0.25))
Change in overlap with self over time
PCA
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 6], ]
indiv_ovlp_l <- lapply(unique(Y$ID), function(x) {
# group data
X <- Y[Y$ID == x, ]
tab_week <- table(X$week)
X <- X[X$week %in% names(tab_week)[tab_week >= 6], ]
if (any(X$week == 1)) {
# run with density
ovlp <- try(rarefact_space_similarity(X = X, dimensions = c("PC1", "PC2"),
group = "week", parallel = 1, pb = FALSE, type = "mean.density.overlap"),
silent = TRUE)
dsts <- try(rarefact_space_similarity(X = X, dimensions = c("PC1", "PC2"),
group = "week", parallel = 1, pb = FALSE, type = "distance"), silent = TRUE)
out <- data.frame(ID = x, group = X$record.group[X$ID == x][1], week = if (is(ovlp,
"try-error"))
NA else ovlp$group.2[ovlp$group.1 == 1], overlap.to.first.week = if (is(ovlp,
"try-error"))
NA else ovlp$mean.overlap[ovlp$group.1 == 1], distance.to.first.week = if (is(dsts,
"try-error"))
NA else dsts$mean.distance[dsts$group.1 == 1], treatment = X$treatment[X$ID ==
x][1], round = X$round[X$ID == x][1])
} else out <- data.frame(ID = x, group = X$record.group[X$ID == x][1], week = NA,
overlap.to.first.week = NA, distance.to.first.week = NA, treatment = X$treatment[X$ID ==
x][1], round = X$round[X$ID == x][1])
return(out)
})
indiv_ovlp <- do.call(rbind, indiv_ovlp_l)
saveRDS(indiv_ovlp, "./data/processed/acoustic_space_density_overlap_to_first_week_by_individual.RDS")indiv_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_first_week_by_individual.RDS")
indiv_ovlp$treatment <- factor(indiv_ovlp$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
indiv_ovlp <- indiv_ovlp[!is.na(indiv_ovlp$group), ]
ggplot(data = indiv_ovlp, aes(x = week, y = overlap.to.first.week, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to acoustic space from first week",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
0.25))
ggplot(data = indiv_ovlp, aes(x = week, y = distance.to.first.week, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Individual distance to acoustic space from first week",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
0.25))
agg_dist <- aggregate(distance.to.first.week ~ week + treatment, indiv_ovlp, mean)
agg_dist$se <- aggregate(distance.to.first.week ~ week + treatment, indiv_ovlp, standard_error)[,
3]
agg_dist$Week <- 1:4
# ggplot(data = agg_dist, aes(x = Week, y = distance.to.first.week, fill =
# treatment)) + geom_bar(stat = 'identity') +
# geom_errorbar(aes(ymin=distance.to.first.week - se, ymax =
# distance.to.first.week + se), width=.1) + # geom_line(size = 1.2) +
# scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~ treatment, ncol =
# 3, scale = 'fixed') + labs(y = 'Distance to self first week', x = 'Week') +
# theme_classic(base_size = 30) + # scale_x_discrete( labels = c(1:4)) +
# theme(legend.position = 'none')
agg_ovlp <- aggregate(overlap.to.first.week ~ week + treatment, indiv_ovlp, mean)
agg_ovlp$se <- aggregate(overlap.to.first.week ~ week + treatment, indiv_ovlp, standard_error)[,
3]
agg_ovlp$Week <- 1:4
# ggplot(data = agg_ovlp, aes(x = Week, y = overlap.to.first.week, fill =
# treatment)) + geom_bar(stat = 'identity') +
# geom_errorbar(aes(ymin=overlap.to.first.week - se, ymax =
# overlap.to.first.week + se), width=.1) + # geom_line(size = 1.2) +
# scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~ treatment, ncol =
# 3, scale = 'fixed') + labs(y = 'Overlap to \nself first week', x = 'Week') +
# theme_classic(base_size = 30) + # scale_x_discrete( labels = c(1:4)) +
# theme(legend.position = 'none')TSNE
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 2], ]
indiv_ovlp_l <- lapply(unique(Y$ID), function(x) {
# group data
X <- Y[Y$ID == x, ]
tab_week <- table(X$week)
X <- X[X$week %in% names(tab_week)[tab_week >= 6], ]
if (any(X$week == 1)) {
# run with density
ovlp <- try(rarefact_space_similarity(X = X, dimensions = c("TSNE1", "TSNE2"),
group = "week", parallel = 1, pb = FALSE, type = "mean.density.overlap"),
silent = TRUE)
dsts <- try(rarefact_space_similarity(X = X, dimensions = c("TSNE1", "TSNE2"),
group = "week", parallel = 1, pb = FALSE, type = "distance"), silent = TRUE)
out <- data.frame(ID = x, group = X$record.group[X$ID == x][1], week = if (is(ovlp,
"try-error"))
NA else ovlp$group.2[ovlp$group.1 == 1], overlap.to.first.week = if (is(ovlp,
"try-error"))
NA else ovlp$mean.overlap[ovlp$group.1 == 1], distance.to.first.week = if (is(dsts,
"try-error"))
NA else dsts$mean.distance[dsts$group.1 == 1], treatment = X$treatment[X$ID ==
x][1], round = X$round[X$ID == x][1])
} else out <- data.frame(ID = x, group = X$record.group[X$ID == x][1], week = NA,
overlap.to.first.week = NA, distance.to.first.week = NA, treatment = X$treatment[X$ID ==
x][1], round = X$round[X$ID == x][1])
return(out)
})
indiv_ovlp <- do.call(rbind, indiv_ovlp_l)
saveRDS(indiv_ovlp, "./data/processed/acoustic_space_density_overlap_to_first_week_by_individual_TSNE.RDS")indiv_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_first_week_by_individual_TSNE.RDS")
indiv_ovlp$treatment <- factor(indiv_ovlp$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
indiv_ovlp <- indiv_ovlp[!is.na(indiv_ovlp$group), ]
ggplot(data = indiv_ovlp, aes(x = week, y = overlap.to.first.week, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to acoustic space from first week",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
0.25))
ggplot(data = indiv_ovlp, aes(x = week, y = distance.to.first.week, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Individual distance to acoustic space from first week",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
0.25))
agg_dist <- aggregate(distance.to.first.week ~ week + treatment, indiv_ovlp, mean)
agg_dist$se <- aggregate(distance.to.first.week ~ week + treatment, indiv_ovlp, standard_error)[,
3]
# agg_stress$Week <- 1:4
# ggplot(data = agg_dist, aes(x = week, y = distance.to.first.week, fill =
# treatment)) + geom_bar(stat = 'identity') +
# geom_errorbar(aes(ymin=distance.to.first.week - se, ymax =
# distance.to.first.week + se), width=.1) + # geom_line(size = 1.2) +
# scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~ treatment, ncol =
# 3, scale = 'fixed') + labs(y = 'Distance to first week', x = 'Week') +
# theme_classic(base_size = 30) + # scale_x_discrete( labels = c(1:4)) +
# theme(legend.position = 'none')Statistical analyses
- on t-SNE dimensions
- predictors were mean-centered
- link functions: binomial distribution for overlap and gaussian for area and call count
- cort.stress and cort.baseline were dropped as they are highly correlated to stress response
- Not a lot of data makes having many predictors complicated
- No random factor (only 1 row per individual)
- No p values printed. Look at whether credibility intervals (l-95% CI & u-95% CI) include 0 to assess statistical significance
Overlap to self (week 1) as response
overlap.to.first.week ~ treatment + stress.response
sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
indiv_ovlp$ID.week <- paste(indiv_ovlp$ID, indiv_ovlp$week, sep = "-")
# addd stress parameters
indiv_ovlp$cort.baseline <- sapply(indiv_ovlp$ID.week, function(x) sels$cort.baseline[sels$ID.week ==
x][1])
indiv_ovlp$cort.stress <- sapply(indiv_ovlp$ID.week, function(x) sels$cort.stress[sels$ID.week ==
x][1])
indiv_ovlp$breath.count <- sapply(indiv_ovlp$ID.week, function(x) sels$breath.count[sels$ID.week ==
x][1])
indiv_ovlp$stress.response <- indiv_ovlp$cort.stress - indiv_ovlp$cort.baseline
# mean centering
indiv_ovlp$cort.baseline <- mean(indiv_ovlp$cort.baseline, na.rm = TRUE) - indiv_ovlp$cort.baseline
indiv_ovlp$cort.stress <- mean(indiv_ovlp$cort.stress, na.rm = TRUE) - indiv_ovlp$cort.stress
indiv_ovlp$breath.count <- mean(indiv_ovlp$breath.count, na.rm = TRUE) - indiv_ovlp$breath.count
indiv_ovlp$stress.response <- mean(indiv_ovlp$stress.response, na.rm = TRUE) - indiv_ovlp$stress.response
# only individuals present in week 1 indiv_ovlp <- indiv_ovlp[indiv_ovlp$ID
# %in% indiv_ovlp$ID[indiv_ovlp$week == 1], ]
# indiv_ovlp <- indiv_ovlp[complete.cases(indiv_ovlp), ]
indiv_ovlp$overlap.to.first.week[indiv_ovlp$overlap.to.first.week == 0] <- 1e-04
indiv_ovlp$overlap.to.first.week[indiv_ovlp$overlap.to.first.week == 1] <- 0.9999week 2:
## week 2
overlap_week_2 <- brm(iter = 5000, silent = 2, overlap.to.first.week ~ treatment +
stress.response, data = indiv_ovlp[indiv_ovlp$week == 2, ], family = Beta(),
control = list(adapt_delta = 0.9), chains = 1)
saveRDS(overlap_week_2, "./data/processed/brms_model_overlap_week_2.RDS")overlap_week_2 <- readRDS("./data/processed/brms_model_overlap_week_2.RDS")
summary_brm_model(overlap_week_2)| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| treatmentMediumStress | 0.065 | 1 | 1527.558 | 5000 | 1 | -0.723 | 0.798 |
| treatmentHighStress | 0.618 | 1 | 1453.076 | 5000 | 1 | -0.131 | 1.334 |
| stress.response | -0.006 | 1 | 2615.287 | 5000 | 1 | -0.026 | 0.013 |
week 4:
## week 4
overlap_week_4 <- brm(iter = 5000, silent = 2, overlap.to.first.week ~ treatment +
stress.response, data = indiv_ovlp[indiv_ovlp$week == 4, ], family = Beta(),
control = list(adapt_delta = 0.9, max_treedepth = 12), chains = 1)
saveRDS(overlap_week_4, "./data/processed/brms_model_overlap_week_4.RDS")overlap_week_4 <- readRDS("./data/processed/brms_model_overlap_week_4.RDS")
summary_brm_model(overlap_week_4)| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| treatmentMediumStress | -0.702 | 1 | 2257.376 | 5000 | 1 | -1.971 | 0.613 |
| treatmentHighStress | 0.454 | 1 | 2095.748 | 5000 | 1 | -0.639 | 1.556 |
| stress.response | 0.042 | 1 | 2357.324 | 5000 | 1 | 0.001 | 0.082 |
Acoustic space change (compare to week 1)
area ~ treatment + stress.response
areas_by_week$ID.week <- paste(areas_by_week$ID, areas_by_week$week, sep = "-")
# addd stress parameters
areas_by_week$cort.baseline <- sapply(areas_by_week$ID.week, function(x) sels$cort.baseline[sels$ID.week ==
x][1])
areas_by_week$cort.stress <- sapply(areas_by_week$ID.week, function(x) sels$cort.stress[sels$ID.week ==
x][1])
areas_by_week$breath.count <- sapply(areas_by_week$ID.week, function(x) sels$breath.count[sels$ID.week ==
x][1])
areas_by_week$stress.response <- areas_by_week$cort.stress - areas_by_week$cort.baseline
# mean centering
areas_by_week$cort.baseline <- mean(areas_by_week$cort.baseline, na.rm = TRUE) -
areas_by_week$cort.baseline
areas_by_week$cort.stress <- mean(areas_by_week$cort.stress, na.rm = TRUE) - areas_by_week$cort.stress
areas_by_week$breath.count <- mean(areas_by_week$breath.count, na.rm = TRUE) - areas_by_week$breath.count
areas_by_week$stress.response <- mean(areas_by_week$stress.response, na.rm = TRUE) -
areas_by_week$stress.response
# only individuals present in week 1 areas_by_week <-
# areas_by_week[areas_by_week$ID %in% areas_by_week$ID[areas_by_week$week ==
# 1], ]
# areas_by_week <- areas_by_week[complete.cases(areas_by_week), ]week 2:
## week 2
area_week_2 <- brm(iter = 5000, silent = 2, size ~ treatment + stress.response, data = areas_by_week[areas_by_week$week ==
2, ], control = list(adapt_delta = 0.9), chains = 1)
saveRDS(area_week_2, "./data/processed/brms_model_area_week_2.RDS")area_week_2 <- readRDS("./data/processed/brms_model_area_week_2.RDS")
summary_brm_model(area_week_2)| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| treatmentMediumStress | 0.175 | 1 | 1787.606 | 5000 | 1 | -0.314 | 0.672 |
| treatmentHighStress | 0.126 | 1 | 1701.261 | 5000 | 1 | -0.384 | 0.614 |
| stress.response | 0.012 | 1 | 2390.428 | 5000 | 1 | -0.002 | 0.026 |
week 4:
## week 4
area_week_4 <- brm(iter = 5000, silent = 2, size ~ treatment + stress.response, data = areas_by_week[areas_by_week$week ==
4, ], control = list(adapt_delta = 0.9), chains = 1)
saveRDS(area_week_4, "./data/processed/brms_model_area_week_4.RDS")area_week_4 <- readRDS("./data/processed/brms_model_area_week_4.RDS")
summary_brm_model(area_week_4)| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| treatmentMediumStress | 0.461 | 1.004 | 1131.119 | 5000 | 1 | -0.335 | 1.255 |
| treatmentHighStress | 0.669 | 1.006 | 1362.654 | 5000 | 1 | -0.054 | 1.394 |
| stress.response | 0.017 | 1.001 | 1806.867 | 5000 | 1 | -0.013 | 0.045 |
Call count
week 2:
## week 2
count_week_2 <- brm(iter = 5000, silent = 2, selec ~ treatment + stress.response,
data = call_count[call_count$week == 2, ], control = list(adapt_delta = 0.9),
chains = 1)
saveRDS(count_week_2, "./data/processed/brms_model_count_week_2.RDS")count_week_2 <- readRDS("./data/processed/brms_model_count_week_2.RDS")
summary_brm_model(count_week_2)| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| treatmentMediumStress | 355.784 | 1.003 | 1074.230 | 5000 | 1 | -127.180 | 804.368 |
| treatmentHighStress | 613.795 | 1.001 | 1050.304 | 5000 | 1 | 131.824 | 1086.057 |
| stress.response | 0.971 | 1 | 1396.041 | 5000 | 1 | -9.032 | 10.761 |
week 4:
## week 4
count_week_4 <- brm(iter = 5000, silent = 2, selec ~ treatment + stress.response,
data = call_count[call_count$week == 4, ], control = list(adapt_delta = 0.9),
chains = 1)
saveRDS(count_week_4, "./data/processed/brms_model_count_week_4.RDS")count_week_4 <- readRDS("./data/processed/brms_model_count_week_4.RDS")
summary_brm_model(count_week_4)| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| treatmentMediumStress | -103.093 | 1 | 1001.866 | 5000 | 1 | -335.787 | 123.953 |
| treatmentHighStress | 3.989 | 1.001 | 840.242 | 5000 | 1 | -352.231 | 359.975 |
| stress.response | 2.278 | 1.001 | 897.594 | 5000 | 1 | -1.867 | 5.900 |
Takeaways
- Stress has an effect on number of calls: high-stress individuals produce more calls than low-stress individuals in week 2
- Stress response has a positive effect on the overlap to week 1 (overlap to itself)
Proposed Models for Stress and Learning Pilot Data based on 9/28/21 discussion
Model set 1: effects of treatment
Predicted variables: - Physiology: Baseline CORT, stress response CORT (stress), difference baseline-stress response CORT (stress response), weight, breath rate - Learning: number of calls, total acoustic area, distance moved from week 1, overlap to self week 1, overlap self to group
General form: - Predicted ~ treatment + week + ind(rand) - Predicted ~ treatment + week + treat * week + ind(rand) - Predicted ~ round + week [to check for effect of round]
Note: Use week 1 as baseline, comparing weeks 2, 3, 4
agg_dat <- aggregate(selec ~ ID + week, data = sels, length)
# compare to week 1 agg_dat$call.count <- sapply(1:nrow(agg_dat), function(x) {
# baseline <- agg_dat$selec[agg_dat$week == 1 & agg_dat$ID == agg_dat$ID[x]] if
# (length(baseline) > 0) change <- agg_dat$selec[x] - baseline else change <-
# agg_dat$selec[x] return(change) } )
# without comparing to week 1
agg_dat$call.count <- sapply(1:nrow(agg_dat), function(x) agg_dat$selec[x])
agg_dat$selec <- NULL
agg_dat$baseline.CORT <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$cort.baseline[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$cort.baseline[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$stress.response <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$stress.response[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$stress.response[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$stress.CORT <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$cort.stress[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$cort.stress[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$weight <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$weight[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$weight[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$breath.rate <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$breath.count[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$breath.count[sels$week == agg_dat$week[x] & sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$acoustic.area <- sapply(1:nrow(agg_dat), function(x) {
area <- areas_by_week$raref.area[areas_by_week$ID == agg_dat$ID[x] & areas_by_week$week ==
agg_dat$week[x]]
if (length(area) < 1)
area <- NA
return(area)
})
agg_dat$acoustic.distance <- sapply(1:nrow(agg_dat), function(x) {
distance <- indiv_ovlp$distance.to.first.week[indiv_ovlp$ID == agg_dat$ID[x] &
indiv_ovlp$week == agg_dat$week[x]]
if (length(distance) < 1)
distance <- NA
return(distance)
})
agg_dat$acoustic.overlap <- sapply(1:nrow(agg_dat), function(x) {
overlap <- indiv_ovlp$overlap.to.first.week[indiv_ovlp$ID == agg_dat$ID[x] &
indiv_ovlp$week == agg_dat$week[x]]
if (length(overlap) < 1)
overlap <- NA
return(overlap)
})
agg_dat$acoustic.overlap.to.group <- sapply(1:nrow(agg_dat), function(x) {
overlap <- group_ovlp$overlap.to.group[group_ovlp$ID == agg_dat$ID[x] & group_ovlp$week ==
agg_dat$week[x]]
if (length(overlap) < 1)
overlap <- NA
return(overlap)
})
agg_dat$treatment <- sapply(1:nrow(agg_dat), function(x) unique(sels$treatment[sels$ID ==
agg_dat$ID[x]]))
agg_dat$round <- sapply(1:nrow(agg_dat), function(x) unique(sels$round[sels$ID ==
agg_dat$ID[x]]))
# remove week 1
agg_dat <- agg_dat[agg_dat$week != 1, ]
responses <- setdiff(names(agg_dat), c("ID", "week", "treatment", "round"))
predictors <- c("~ treatment + week + (1|ID)", "~ treatment * week + (1|ID)", "~ round + week + (1|ID)")
formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)vars_to_scale <- names(agg_dat)[sapply(agg_dat, is.numeric)]
vars_to_scale <- vars_to_scale[!vars_to_scale %in% c("week", "round")]
for (i in vars_to_scale) agg_dat[, vars_to_scale] <- scale(agg_dat[, vars_to_scale])
treatment_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
]
sub_dat
mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
chains = 2), silent = TRUE)
return(mod)
})
names(treatment_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(treatment_models, "./data/processed/treatment_models.RDS")
agg_dat$week <- as.factor(agg_dat$week)
treatment_models_week_factor <- lapply(1:nrow(formulas), function(x) {
sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
]
mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
chains = 2), silent = TRUE)
return(mod)
})
names(treatment_models_week_factor) <- paste(formulas$responses, formulas$predictors)
saveRDS(treatment_models_week_factor, "./data/processed/treatment_models_week_factor.RDS")week as continuous
treatment_models <- readRDS("./data/processed/treatment_models.RDS")
for (x in treatment_models) {
summary_brm_model(x)
}| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | -0.342 | 1.001 | 2032.331 | 5000 | 2 | -0.623 | -0.047 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | 0.057 | 1 | 3386.12 | 5000 | 2 | -0.198 | 0.305 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | -0.08 | 1.001 | 2172.071 | 5000 | 2 | -0.368 | 0.209 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | 0.018 | 1.001 | 3184.54 | 5000 | 2 | -0.28 | 0.314 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | 0.062 | 1 | 2669.238 | 5000 | 2 | -0.227 | 0.366 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | 0.267 | 1.001 | 1967.507 | 5000 | 2 | -0.04 | 0.569 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | -0.14 | 1 | 3307.176 | 5000 | 2 | -0.397 | 0.111 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | 0.067 | 1 | 2819.069 | 5000 | 2 | -0.213 | 0.339 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | -0.089 | 1 | 2859.569 | 5000 | 2 | -0.384 | 0.212 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | 0.001 | 1 | 3164.738 | 5000 | 2 | -0.309 | 0.302 |
Takeaways
- No significant nor marginal effects
Model set 2: individual stress and effects on learning
General form: - Learning variables ~ baseline CORT [Week 2 immediate stress and learning, week 3 chronic stress and learn]
agg_dat <- agg_dat[agg_dat$week != "4", ]
agg_dat$week <- as.factor(agg_dat$week)
responses <- c("call.count", "acoustic.area", "acoustic.distance", "acoustic.overlap",
"acoustic.overlap.to.group")
predictors <- c("~ baseline.CORT + week + (1|ID)", "~ stress.response + week + (1|ID)")
formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)treatment_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
]
mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
chains = 1), silent = TRUE)
return(mod)
})
names(treatment_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(treatment_models, "./data/processed/stress_on_learning_models.RDS")treatment_models <- readRDS("./data/processed/stress_on_learning_models.RDS")
for (x in treatment_models) {
summary_brm_model(x)
}| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | -0.192 | 1.002 | 997.111 | 5000 | 1 | -0.546 | 0.177 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | -0.154 | 1 | 537.359 | 5000 | 1 | -0.462 | 0.159 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | 0.007 | 1.015 | 256.643 | 5000 | 1 | -0.248 | 0.274 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | 0.101 | 1 | 616.351 | 5000 | 1 | -0.247 | 0.44 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | 0.299 | 1.001 | 1394.371 | 5000 | 1 | -0.112 | 0.69 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | 0.013 | 1 | 556.089 | 5000 | 1 | -0.25 | 0.287 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | -0.058 | 1.001 | 677.164 | 5000 | 1 | -0.391 | 0.286 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | -0.138 | 1.002 | 869.593 | 5000 | 1 | -0.49 | 0.209 |
Takeaways
- No significant effects
Model set 3: individual stress and effects on learning on week 2 and 3 independently
Week 2 immediate stress and learning, week 3 chronic stress and learn
General form: - Learning variables ~ baseline CORT - Learning variables ~ difference baseline-stress CORT
Week 2
agg_dat_w2 <- agg_dat[agg_dat$week == "2", ]
responses <- c("call.count", "acoustic.area", "acoustic.distance", "acoustic.overlap",
"acoustic.overlap.to.group")
predictors <- c("~ baseline.CORT + (1|ID)", "~ stress.response + (1|ID)")
formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)
treatment_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- agg_dat_w2[!is.na(agg_dat_w2[names(agg_dat_w2) == formulas$responses[x]]),
]
mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
chains = 1), silent = TRUE)
return(mod)
})
names(treatment_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(treatment_models, "./data/processed/stress_on_learning_models_week_2.RDS")treatment_models <- readRDS("./data/processed/stress_on_learning_models_week_2.RDS")
for (x in treatment_models) {
summary_brm_model(x)
}| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | -0.177 | 1 | 1201.316 | 5000 | 1 | -0.501 | 0.164 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | 0.037 | 1.002 | 701.318 | 5000 | 1 | -0.328 | 0.408 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | -0.138 | 1.012 | 135.412 | 5000 | 1 | -0.426 | 0.133 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | -0.03 | 1.045 | 78.7 | 5000 | 1 | -0.329 | 0.298 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | 0.11 | 1.004 | 1054.756 | 5000 | 1 | -0.257 | 0.458 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | 0.261 | 1.001 | 831.794 | 5000 | 1 | -0.072 | 0.586 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | -0.178 | 1 | 227.403 | 5000 | 1 | -0.545 | 0.202 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | 0.013 | 1 | 1041.844 | 5000 | 1 | -0.232 | 0.25 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | -0.045 | 1.002 | 546.832 | 5000 | 1 | -0.381 | 0.267 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | -0.138 | 1.004 | 1319.857 | 5000 | 1 | -0.529 | 0.261 |
Week 3
agg_dat_w3 <- agg_dat[agg_dat$week == "3", ]
responses <- c("call.count", "acoustic.area", "acoustic.distance", "acoustic.overlap",
"acoustic.overlap.to.group")
predictors <- c("~ baseline.CORT + (1|ID)", "~ stress.response + (1|ID)")
formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)
treatment_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- agg_dat_w3[!is.na(agg_dat_w3[names(agg_dat_w3) == formulas$responses[x]]),
]
mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
chains = 1), silent = TRUE)
return(mod)
})
names(treatment_models) <- paste(formulas$responses, formulas$predictors)
saveRDS(treatment_models, "./data/processed/stress_on_learning_models_week_3.RDS")treatment_models <- readRDS("./data/processed/stress_on_learning_models_week_3.RDS")
for (x in treatment_models) {
summary_brm_model(x)
}| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | -0.192 | 1.002 | 997.111 | 5000 | 1 | -0.546 | 0.177 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | -0.154 | 1 | 537.359 | 5000 | 1 | -0.462 | 0.159 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | 0.007 | 1.015 | 256.643 | 5000 | 1 | -0.248 | 0.274 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| baseline.CORT | 0.101 | 1 | 616.351 | 5000 | 1 | -0.247 | 0.44 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | 0.299 | 1.001 | 1394.371 | 5000 | 1 | -0.112 | 0.69 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | 0.013 | 1 | 556.089 | 5000 | 1 | -0.25 | 0.287 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | -0.058 | 1.001 | 677.164 | 5000 | 1 | -0.391 | 0.286 |
| Estimate | Rhat | Bulk_ESS | iterations | chains | CI_low | CI_high | |
|---|---|---|---|---|---|---|---|
| stress.response | -0.138 | 1.002 | 869.593 | 5000 | 1 | -0.49 | 0.209 |
Takeaways
- No significant nor marginal effects
Session information
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=pt_BR.UTF-8
## [7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] PhenotypeSpace_0.1.0 warbleR_1.1.26 NatureSounds_1.0.4
## [4] seewave_2.1.8 tuneR_1.3.3.1 ggridges_0.5.3
## [7] posterior_1.2.0 ggrepel_0.9.1 cowplot_1.1.1
## [10] tidybayes_3.0.2 modelr_0.1.8 tidyr_1.2.0
## [13] forcats_0.5.1 purrr_0.3.4 dplyr_1.0.8
## [16] lme4_1.1-28 Matrix_1.3-4 brms_2.16.3
## [19] Rcpp_1.0.8 vegan_2.5-7 lattice_0.20-44
## [22] permute_0.9-7 raster_3.5-15 spatstat_2.3-3
## [25] spatstat.linnet_2.3-2 spatstat.core_2.4-0 rpart_4.1-15
## [28] nlme_3.1-152 spatstat.random_2.1-0 spatstat.geom_2.3-2
## [31] spatstat.data_2.1-2 GGally_2.1.2 adehabitatHR_0.4.19
## [34] adehabitatLT_0.3.25 CircStats_0.2-6 boot_1.3-28
## [37] adehabitatMA_0.3.14 ade4_1.7-18 deldir_1.0-6
## [40] sp_1.4-6 MASS_7.3-54 Rtsne_0.15
## [43] randomForest_4.7-1 formatR_1.11 knitr_1.37
## [46] kableExtra_1.3.4 ggplot2_3.3.5 viridis_0.6.2
## [49] viridisLite_0.4.0 pbapply_1.5-0 readxl_1.3.1
## [52] lubridate_1.8.0 remotes_2.4.2
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 tidyselect_1.1.2 fftw_1.0-6.1
## [4] htmlwidgets_1.5.4 grid_4.1.0 munsell_0.5.0
## [7] codetools_0.2-18 DT_0.20 miniUI_0.1.1.1
## [10] withr_2.4.3 Brobdingnag_1.2-7 colorspace_2.0-3
## [13] highr_0.9 rstudioapi_0.13 stats4_4.1.0
## [16] dtw_1.22-3 tensor_1.5 bayesplot_1.8.1
## [19] labeling_0.4.2 rstan_2.21.3 polyclip_1.10-0
## [22] farver_2.1.0 bridgesampling_1.1-2 coda_0.19-4
## [25] vctrs_0.3.8 generics_0.1.2 xfun_0.29
## [28] R6_2.5.1 markdown_1.1 HDInterval_0.2.2
## [31] gamm4_0.2-6 projpred_2.0.2 bitops_1.0-7
## [34] spatstat.utils_2.3-0 reshape_0.8.8 assertthat_0.2.1
## [37] promises_1.2.0.1 scales_1.1.1 gtable_0.3.0
## [40] processx_3.5.2 goftest_1.2-3 rlang_1.0.1
## [43] systemfonts_1.0.4 splines_4.1.0 broom_0.7.12
## [46] checkmate_2.0.0 inline_0.3.19 yaml_2.3.5
## [49] reshape2_1.4.4 abind_1.4-5 threejs_0.3.3
## [52] crosstalk_1.2.0 backports_1.4.1 httpuv_1.6.5
## [55] rsconnect_0.8.25 tensorA_0.36.2 tools_4.1.0
## [58] ellipsis_0.3.2 jquerylib_0.1.4 RColorBrewer_1.1-2
## [61] proxy_0.4-26 plyr_1.8.6 base64enc_0.1-3
## [64] RCurl_1.98-1.6 ps_1.6.0 prettyunits_1.1.1
## [67] zoo_1.8-9 cluster_2.1.2 magrittr_2.0.2
## [70] ggdist_3.1.0 colourpicker_1.1.1 mvtnorm_1.1-3
## [73] matrixStats_0.61.0 shinyjs_2.1.0 mime_0.12
## [76] evaluate_0.15 arrayhelpers_1.1-0 xtable_1.8-4
## [79] shinystan_2.5.0 gridExtra_2.3 rstantools_2.1.1
## [82] compiler_4.1.0 tibble_3.1.6 crayon_1.5.0
## [85] minqa_1.2.4 StanHeaders_2.21.0-7 htmltools_0.5.2
## [88] mgcv_1.8-36 later_1.3.0 RcppParallel_5.1.5
## [91] DBI_1.1.1 cli_3.2.0 parallel_4.1.0
## [94] igraph_1.2.11 pkgconfig_2.0.3 signal_0.7-7
## [97] terra_1.5-21 spatstat.sparse_2.1-0 xml2_1.3.3
## [100] svUnit_1.0.6 dygraphs_1.1.1.6 svglite_2.1.0
## [103] bslib_0.3.1 webshot_0.5.2 rvest_1.0.2
## [106] stringr_1.4.0 distributional_0.3.0 callr_3.7.0
## [109] digest_0.6.29 rmarkdown_2.11 cellranger_1.1.0
## [112] shiny_1.7.1 gtools_3.9.2 rjson_0.2.21
## [115] nloptr_2.0.0 lifecycle_1.0.1 jsonlite_1.8.0
## [118] fansi_1.0.2 pillar_1.7.0 loo_2.4.1
## [121] fastmap_1.1.0 httr_1.4.2 pkgbuild_1.3.1
## [124] glue_1.6.1 xts_0.12.1 shinythemes_1.2.0
## [127] stringi_1.7.6 sass_0.4.0